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First-principles accurate total-energy surfaces for polar structural distortions of 
BaTiOs, PbTiOa, and SrTiOs: consequences to structural transition temperatures 

Takeshi Nishimatsu^ , Masaya Iwamoto^, Yoshiyuki Kawazoe^, and Umesh V. Waghmare^ 

^Institute for Materials Research (IMR), Tohoku University, Sendai 980-8577, Japan 

^ Theoretical Sciences Unit, Jawaharlal Nehru Centre for Advanced 

Scientific Research (JNCASR), Jakkur, Bangalore, 560 064, India 

Specific forms of tiie excliange correlation energy functionals in first-principles density functional 
theory-based calculations, such as the local density approximation (LDA) and generalized-gradient 
approximations (GGA), give rise to structural lattice parameters with typical errors of —2 % and 2 %. 
Due to a strong coupling between structure and polarization, the order parameter of ferroelectric 
transitions, they result in large errors in estimation of temperature dependent ferroelectric structural 
transition properties. Here, we employ a recently developed GGA functional of Wu and Cohen 
[Phys. Rev. B, 73, 235116 (2006)] and determine total-energy surfaces for zone-center distortions 
of BaTiOa, PbTiOa, and SrTiOa, and compare them with the ones obtained with calculations 
based on standard LDA and GGA. Confirming that the Wu and Cohen functional allows better 
estimation of structural properties at K, we determine a new set of parameters defining the effective 
Hamiltonian for ferroelectric transition in BaTiOa. Using the new set of parameters, we perform 
molecular-dynamics (MD) simulations under effective pressures p = 0.0 GPa, p = —2.0 GPa, andp — 
— O.OOST GPa. The simulations under p = — O.OOST GPa, which is for simulating thermal expansion, 
show a clear improvement in the cubic to tetragonal transition temperature and c/a parameter 
of its ferroelectric tetragonal phase, while the description of transitions at lower temperatures to 
orthorhombic and rhombohedral phases is marginally improved. Our findings augur well for use of 
Wu-Cohen functional in studies of ferroelectrics at nano-scale, particularly in the form of epitaxial 
films where the properties depend crucially on the lattice mismatch. 
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I. INTRODUCTION 

It is well known that first-principles density functional 
theory based calculations within the local density approx- 
imation (LDA) underestimate lattice constants slightly 
(1-2%), and consequently calculated double-well total- 
energy surfaces-'^ for ferroelectric structural distortions 
of ABO3 perovskite-type ferroelectrics are shallower giv- 
ing the theoretical transitions temperatures much lower 
than their observed values (Tc = 403 K = 0.0347 eV 
for BaTiOa, for example)^. In Monte-Carlo (MC) 
simulations^ and molecular-dynamics (MD) simulations^ 
of BaTiOa, a perovskite-type ferroelectric, an effective 
Hamiltoniar^ii^ constructed from LDA calculations was 
used. To overcome the limitation of underestimation of 
lattice constant, these simulations were carried out with 
a negative pressure of —5 GPa. Similarly, Monte Carlo 
simulations for PbTiOa were carried out with an effective 
Hamiltonian^ constructed from LDA calculations at the 
experimental lattice constant. However, both of these 
schemes resulted in underestimation of Tq ■ 

There are two main sources of errors in estimation of 
Tc in such simulations: (a) neglect of anharmonic cou- 
pling between soft modes and higher energy phonons in 
construction of effective Hamiltonian, and (b) those aris- 
ing from underestimation of lattice constants in DFT cal- 
culations. While the former was assessed to be small 
in earlier work^ and can be partly corrected using a T- 
dependent pressure to yield a correct thermal expansion, 
a systematic investigation of the latter would be useful 



in planning and evaluating future first-principles simula- 
tions of ferroelectrics. 

To overcome the limitation of DFT calculations in es- 
timation of structural parameters in ferroelectrics, Wu 
and Cohen introduced a new flavor of generalized gra- 
dient approximation (GGA), and obtained an excellent 
agreement between calculated and experimentally ob- 
served lattice constants at zero Kelvin for FbTiOa and 
BaTiOa®. Recent theoretical works^iiS further strength- 
ened that the Wu and Cohen functional gives acceptable 
structural properties of ABO3 ferroelectrics such as lat- 
tice constants, c/a ratio, atomic displacements, phonon 
frequencies. Born effective charges, etc. 

Here, we use the Wu and Cohen GGA-functional and 
determine possibly more realistic total-energy surfaces 
of polar distortions of BaTiOg, FbTiOg, and SrTiOg 
than those from LDA calculations. Further, we con- 
struct an effective Hamiltonian for BaTiOs with a set of 
parameters determined from first-principles calculations 
based on Wu-Cohen functional, and estimate the three of 
its transition temperatures, and temperature dependent 
structural properties. Through comparison with transi- 
tions properties obtained with the LDA-based effective 
Hamiltonian and from experiment, we evaluate the effi- 
cacy of Wu-Cohen functional in determination of finite- 
temperature properties. 

In Sec. ini we describe the formalism of methods used 
in computations, with a focus on details of the proce- 
dure for determination of the set of parameters of the 
effective Hamiltonian, which is slightly different from the 



earlier works^. In Sec. IIIIl we present a comparative 
analysis of calculated total-energy surfaces of BaTiOa, 
PbTiOa, and SrTiOa and include results of MD simu- 
lations of BaTiOs, and finally summarize our work and 
conclusions in Sec. HVl 



II. 



METHODS OF CALCULATION AND 
FORMALISM 

A. First-principles methods 



All calculations are performed with ABINIT codeiii^. 
Bloch wave functions of electrons are expanded in terms 
of plane waves with a cut-off energy of 60 Hartree, and 
are sampled on an 8x8x8 grid of fc-points in the first Bril- 
louin zone. We use different choices of exchange correla- 
tion energy functionals. For LDA calculations, we use the 
one parametrized by Teter— along with Teter's extended 
norm-conserving pseudopotentialsi^. For GGA calcu- 
lations, we use "PBE"J^ and "Wu and Cohen"® func- 
tionals, along with Rappe's optimized pseudopotentials^^ 
generated with Opium code-- and compare their results. 



B. Total-energy surface 
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FIG. 1. (Color online) Total-energy surfaces for zone-center 
distortions of a ferroelectric ABO3 perovskite on a two- 
dimensional subspace (u^ ,vf) of the atomic-displacement 
space («2 , «^ , Hz ^ , «z " , «z '" ) . (a) Schematic contour plot for 
atomic displacements from the centrosymmetric cubic struc- 
ture (b), ezz = 0, is compared to (c) that from the tetragonal 
structure (d), Ezz > 0. Thick solid lines are the valley lines for 
fixed e^z's. Dashed lines show the direction of the Fis soft- 
mode eigenvector ^ at zero strain. Note that the direction is 
tangential to the valley line at d^ = for Czz = 0, but this is 
not the case for ezz 7^ 0. 



In 1994, King-Smith and Vanderbilt studied the total- 
energy surface for zone-center distortions of perovskite- 
type ferroelectric oxides ABO3 at zero tempera- 
ture using first-principles calculations with ultrasoft- 
pseudopotentials and a plane-wave basis setji Starting 
from the centrosymmetric cubic perovskite structure, 
and using the normalized Fis soft-mode eigenvector ^^ 
(= ^x = ^y = ^zj due to the cubic symmetry) of the 
interatomic force constant (IFC) matrijii^, they define 
displacements v"^ of atoms t (=A, B, Oi, Oh, Om) in 
the Cartesian directions a(= x,y,z) as 



/ v^ \ 



-^QSa '^ 



V <"' / 



Sq 

Set 
;^Oii 



(1) 



with the scalar soft-mode amplitude Uq. Under the con- 
dition that the strain components r]i (i — I, ■ ■ ■ , 6; Voigt 
notation; 771 = Cxx, Vi — ^yz) minimize the total energy 
for each u = {ux, Uy,Uz), they expressed the total energy 
as 



E' 
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where u^ — Ux + Uy + u1, E^ is the total energy of the cu- 
bic structure, k is half the eigenvalue of the Fis soft mode, 
and a' and 7' are anharmonic coefficients including the 
contribution of relaxation of strain through its coupling 
with atomic displacements. Similar analysis is used in 



constructions of effective Hamiltonian in Refs. Q and m, 
which assumes that anharmonic coupling between soft 
modes and other IR-active modes is vanishingly small. 

In 2004, Hashimoto, Nishimatsu et al. included some 
of these anharmonic effects by redefining Ua as 



ua - y«)' + {v^f + (z;°^)' + {v2"y + [v^'" 

(3) 
and developed automatic computational method to de- 
termine valley line of the total-energy surface in the 15- 
dimensional coordinate space of atomic displacements^ 
(See Fig.[T]). This method reveals that normalized "direc- 
tion" ^ of the atomic displacements from the centrosym- 
metric cubic phase to the distorted minima is not con- 
stant as in Eq. ([2]), but a function of it, i.e. ${u), and the 
total-energy surface generally cannot be expressed with 
a 4th order function of atomic displacement amplitude 
u as in Eq. ([2]). ^(m) is determined by minimizing en- 
ergy with respect to all atomic displacements {w^} such 
that {ua} in Eq. (3) is fixed. In this paper, we em- 
ploy this valley line tracing method as implemented in 
the Patched ABINIT version 5.7.3 to determine accurate 
total-energy surfaces. The patch file for ABINIT is in 
the EPAPSi^. The valley lines can be also calculated 
under any positive or negative pressure, through use of 
enthalpy H = E + pV and correspondingly the enthalpy 
differences H — H'^ should be compared rather than total 
energy E — E^ . 



C. Effective Hamiltonian 

The effective Hamiltonian constructed from first- 
principles calculations and used in MD simulations is ba- 
sically the same as that in Ref. [a and 0, 



R,a R,a 



wliR) 



+ F^°"({m}) + V'^P\{u}) + V'^°'\{u}) 

ycoup,homo(|^|^ r/i, • • •, 776) + V^™"P' '"''°({m}, {w}) 

-Z*Y,£-u{R). (4) 



H 



Detailed explanation of symbols in the effective Hamilto- 
nian can be found in Refs. [l|, 0, and [^- We newly intro- 
duced 6th-order terms with coefficients fci , ^2 , and k^ and 
an 8th-order term k4U^{R) to the local- mode self-energy 
F"°"({m}) as 



R 

+ 7 [ul{R)uliR) + ul{R)ui{R) + ul{R)ul{R)\ 

+ kiu%R) + fc2 [utiR)iul{R) + ul{R)) 

+ u^y{R){ul{R) + ul{R)) + ut{R){ul{R) + ul{R))] 

+ kzul{R)ul{R)ul{R) + fc4w'(fi)} , (5) 

where u'^{R) = ul{R) + ul{R) + ul{R). We introduced 
this 6th-order terms to follow-up the total-energy surface 
precisely. The 8th-order term is for preventing the |m| — > 
oo breakdown under negative fci. 

In next sections IHDI and IHEl we explain how to de- 
termine the parameters for the effective Hamiltonian of 
Eq. (IH) in detail. 

D. Elastic coefficients and total-energy surface 

Elastic constants expressed in energy unit Bn = OqCh 
and Bi2 = a!^Ci2., where ao is the equilibrium lattice con- 
stant in cubic structure, can be calculated by deforming 
the cubic unit cell of ABO^ with strain tensors 



and 



^ 



^ 




(6) 



(7) 



Deformation in Eq. © alters the total energy from its 
equilibrium value £'" by 

E{5)^E'' + ^{Bii+2Bi2)5^ + 0{5^) . (8) 



More precisely, volume dependence of total energy 
may be fitted with the Birch-Murnaghan equation of 
stated"—. Deformation in Eq. ^ gives 



E{5) ^ E"" + -Bii5^ + 0{5^) 



For i?44 = 006*44, deformation 




and 



E{5) 



T^B^' 



0{5^) 



(9) 



(10) 



(11) 



can be used. 

BixxiBiyy, and B^yz, the coupling coefficients defined 
in Ref. 'll, are determined from quadratic u dependence 
of strain. In the case of [110] distortion (see Fig.[5Je), for 
example). 



^xx — ^xx^ 



^xy — ^xy ^ 
^zz — ^zz^ 



emerge 



Bixx — —^Biiaxx + 2(i?ii — 2^12)022 
-^lyy ~ —^Bi2axx — ^BnUzz 

B^yz ^ ~2B/^/^axy ■ 



(12a) 
(12b) 
(12c) 

(13a) 
(13b) 
(13c) 



Anharmonic coefficients in the on-site energy a, 7, 
ki, k2, fcs, and fc4 in Eq. ([5]) are determined from u- 
dependcnces of total energies of [001], [110], and [111] 
distortions as 



Eooi{u) =Ku^ + a'u + kiu^ -t- /c4M^, 

Eiiq{u) =ku^ + [a' + jl')u'^ + (fci + jk2)u^ + k^u^, 



(14a) 
(14b) 



Eiiiiu) =KU^ + {a' + -Y)u* 
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+ (fcl + gfc2 + ^^3)"*^ + hU^ 



(14c) 



With Eq. (19a) and (19b) in Ref. [ij, a' and 7' can be 
converted into a and 7. It should be mentioned that 
it is quite difficult to express the total-energy surfaces 
even with up to 8th order polynomial in wide range of u. 
Therefore, we fit Eqs. (J14al) - (|14cp only to the calculated 
data points within narrow range of u, e.g. |u| < 0.3 [A] 
for BaTiOs. 



E. Response-function calculations 

We perform some response-function (RE) 
calculations^'^ with ABINIT, determine IFC matri- 
ces at the fc-points of F, X, M, R, and center of the S 



axis (See Fig.[3lA)), then calculate their eigenvalues and 
eigenvectors. 

We can determine local and short-range interaction 
parameters K2 and ji , • • • , J7 in Ref. |j from selected 
eigenvalues 2K(rTo), 2k(Xi), 2k(X5), 2k(M3-), 2k(M5.), 

I 



2k(R25'), and 2k(Elo)- Here, it is emphasised that K{ki) 
is half of the mode-z eigenvalue 2«;(fci) of the IFC matrix 
at each fc-point. Practically, K2 and ji, ■ • • , jV are deter- 
mined by solving linear equation as described in Ref. y, 
in CGS units, 
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(15b) 

(15c) 

(15d) 

(15e) 
(15f) 
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and 



= J6 - J7 



(16) 



Moreover, we newly assume that J5 = and jj — 0. With 
this assumption, we can omit one RF calculation for the 
center of the S axis, and we do not use Eq. ( |15gP and 
Eq. ([T6|) . We can employ this assumption because J5 and 
J7 do not affect low energy polarization modes. 

There may be an inconsistency between k from 
Eqs. (fT4a|) - ((T4c)) and K(rTo) in Eq. (fT5a| . To keep the 
total-energy surfaces unchanged, k from Eqs. (|14ap - (|14cp 
should be adopted. Therefore, We add the difference be- 
tween K and /{(Fto) to K2 as 



K2 ^ K2 + [k- k(Fto)] 



(17) 



This correction can be employed, because correction of K2 
just leads parallel elevation of dispersion, e.g. Fig. [31JB). 
From the calculated normalized Fis soft-mode eigen- 
vector ^a , we determine the Born effective charge Z* and 
the effective mass M^jp^j^, of the soft mode. For Z* , we 
also use the calculated effective charge tensor Z*^ for 
each atom r, 



r 

The effective mass will be 



(18) 



(19) 



where M"^ is the mass of atom t. 

The optical dielectric constant e^o can be also deter- 
mined in RF calculations. 



F. Conditions of molecular-dynamics simulations 

MD simulations for BaTi03 with the effective Hamilto- 
nian of Eq. (|4]) are performed with our original MD code 
fercun (http://loto.sourceforge.net/fercmi/). De- 
tails of the code can be found in Ref. 5. Temperature 
is kept constant in each temperature step in the canoni- 
cal ensemble using the Nose-Poincare thermostat. ^"^ This 
simplectic thermostat is so efficient that we can set the 
time step to At == 2 fs. In our present MD simulations, we 
thermalize the system for 180,000 time steps, after which 
we average the properties for 20,000 time steps. We used 
a supercell of system size L^, x Lj, x L^ = 14 x 14 x 14 
and small temperature steps in heating-up (-1-1 K/step) 
and cooling-down (—1 K/step) simulations. It should be 
noted that the larger supercell size and the more rapid 
heating-up and cooling-down result in the larger temper- 
ature hysteresis. The initial configuration are generated 
randomly: (u^) = 0.11 A and {uD - {u^}'^ = (0.02 A)^. 
We have checked that there is no dependence of results 
of these simulations on initial configurations. 



TABLE I. Calculated equilibrium lattice constants ao for cu- 
bic phases are compared to experimental values observed just 
above Tc- Experimental values are cited from Refs. [2J andQ 
for BaTiOa (Tc = 403 K), Ref.lH for PbTiOs (Tc = 763 K), 
and Ref.lil for SrTiOs (Tc = 106 K). 



III. RESULTS AND DISCUSSION 

Calculated total-energy curves along [001], [110], and 
[111] distortion directions of BaTiOs with three function- 
als are shown in Fig. [2] It is clear that LDA results in 
shallow double wells and GGA (PBE) results in rather 
deep wells, while GGA (Wu and Cohen) results in in- 
termediate depths of the double well potentials. Note 
that, in Fig. EJa), double wells of LDA results cannot be 
recognized in this scale of energy. It can be said that 
GGA (Wu and Cohen) succeeds in reproducing total- 
energy surfaces at K. These results can be understood 
from estimated equilibrium cubic lattice constants Oq — 
3.938, 4.034, and 3.986 A, respectively (See also TableU). 
GGA (Wu and Cohen) results also reveal that the "di- 
rection" ^(m) of atomic displacements largely depends on 
u, as shown in Fig. [2jf) even in BaTiOa, which exhibits 
relatively small polar structural distortions across its fer- 
roelectric transition. For emphasizing that the depth of 
the double wells are strongly affected by equilibrium cu- 
bic lattice constant, total-energy surfaces calculated with 
LDA under negative pressure —7.0 GPa (oq — 3.989 A) 
are calculated and show in Fig. ^d) . This value of neg- 
ative pressure, —7.0 GPa, is selected to get similar total- 
energy surfaces as the GGA (Wu and Cohen) results from 
— 1.0, — 2.0, • ■ • ,—9.0 GPa calculations. It can be seen 
that LDA under certain negative pressure gives results 
similar to those of GGA (Wu and Cohen). 

From our calculations with GGA (Wu and Cohen), i.e. 
total-energy surfaces, u-dependence of strain (Fig.[2Ke)), 
IFC matrices, etc., we construct a new parameter set 
of effective Hamiltonian for BaTiOa. The parameters 
from Refs. Ill and 4 and those of present work are com- 
pared in Table HIl As shown in Fig.[2IA), without short- 
range interaction, pure dipole-dipole long-range interac- 
tion results in an antiferroelectric cell-doubling state as 
the most stable structure, corresponding to the strongest 
instability at M point. However, as shown in Fig. |3KB), 
introduction of short-range interactions K2 and j'l , • • • , J7 
results in the ferroelectric state as the most stable struc- 
ture at the r point. 

In Fig. 21 dipole moment per unit cell as a function of 
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FIG. 2. (Color online) (a)-(d) Total-energy surfaces for zone- 
center distortions of BaTiOs . (e) u-dependence of strain along 
[110] distortion, (f) "Direction" ^(n) of atomic displacements 
along [001] distortion. 




TABLE II. Comparison of two set of parameters for the 
BaTiOs effective Hamiltonian. " — " indicates that values was 
not in use. "n.a." indicates that values are not available, p 
is constant or temperature T [K] dependent effective negative 
pressures appUed while MD simulations, k in Eqs. (|14ap - 
(|14c[) and niki) in Eqs. H15ap -( [T5g| are also listed. They are 
used to determine K2 and ji,--- ,J7. 51 are the soft-mode 
eigenvector. 
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FIG. 3. (Color onhne) (A) Half of eigenvalues of the 3 x 
3 long-range dipole-dipole interaction matrix $(fc) (Fourier 
transform of Eq. (10) in Ref . [g) are plotted along symmetric 
axes in the the first Brillouin zone of the simple-cubic lattice. 
Special points and fc/(27r/a) = (i, |,0) (the center of the E 
axis) are indicated with vertical dotted lines. Labels (a)-(g) 
corresponds to Eqs. (|15a|l -( [T5g| , respectively. Tics in the unit 
of — — J is placed in left side. Tics in the unit of eV, in the 

case of the parameter set of Table |IIJ is placed in right side. 
(B) Half of eigenvalues of the total (long-range -|- short-range) 
interaction matrix $'>"^'^(fc) (Eq. (13) in Ref.|H). 



u for atomic displacements along [001] distortion calcu- 
lated with the Berry-phase theory^^ is shown. As com- 
pared with Z*M, it can be seen that linearity is broken 
at large u. Since following MD simulations are based 
in energetics, effects of this nonlinear Z*(u) mainly get 
folded into anharmonic terms, i.e. a, 7, fci, ^2, fcs, and 
^4. However, still, there might be issues with intersite 
anharmonic interactions. We leave the issues for future 
studies. 

From the Table HIl it is evident that most parameters 
in the effective Hamiltonian are sensitive to the choice 
of exchange-correlation functional, with the exception of 
elastic constants and the mode effective charge. Largest 
change is seen in the parameters of coupling between 
strain and polarization, as is expected from the fact that 
Wu-Cohen functional gives a better estimate of lattice 
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FIG. 4. (Color online) Calculated dipole moment per unit 
cell as a function of u for atomic displacements along [001] 
distortion (solid line). Z*u is also plotted for comparison 
(dashed line). 



parameters and that they couple strongly with polariza- 
tion. Electronic contribution to the dielectric constant 
too is further overestimated with Wu-Cohen functional 
(it is typically 20 % overestimated in an LDA-based cal- 
culation) . Due to the use of Hashimoto-Nishimatsu's val- 
ley tracing technique, description of the on-site potential 
energy curve requires anharmonic terms expanded up to 
8th order. 

Using the Wu-Cohen functional-based parametrized ef- 
fective Hamiltonian, we perform heating-up and cooling- 
down MD simulations. In Fig. [S] lattice parameters 
as functions of temperature are plotted under (a) p = 
0.0 GPa, (b) p = -2.0 GPa, and (c) v = -0.005T GPa, 
where T is the temperature in Kelvin, p — —2.0 GPa 
and p = — 0.005r are used to obtain a lattice con- 
stant of 4.010 A just above Tq = 403 K. Note that, 
at T = 400 K, p = -O.OOST ^ -2.0 GPa. The lat- 
ter temperature-dependent effective negative pressure is 
for simulating thermal expansiou i^^'^^ We note that ad- 
justment of the lattice constant through a negative pres- 
sure(s) results in significant improvement in the temper- 
ature of transition from cubic to tetragonal phase, while 
the description of the lower two transitions improves only 
slightly. Results with Wu-Cohen functional (in all the 
three schemes, 0.0, —2.0, and — O.OOST GPa pressures) 
show a significant improvements in cubic-to-tetragonal 
transition temperature compared to the LDA-based re- 
sults of previous MC^ and MD'^ calculations, but not in 
the lower two transitions. These mal-improvements in 
the lower two transitions may be coming from difficulties 
in accurate first-principles calculations and polynomial 
fittings of almost degenerated bottom of double wells of 
[001], [110], and [111] distortions. InTablelTTIl simulated 
transition temperatures are compared to the previous 
MD simulations'^ with LDA-based parameters and exper- 
imentally observed values. Comparing them to experi- 
mentally observed temperature dependence of lattice pa- 
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FIG. 5. (Color online) Simulated temperature dependence 
of lattice parameters under (a) p = —0.0 GPa, (b) p = 
-2.0 GPa, (c) V = -O.OOST GPa. 



rameters in Ref. JJ, we also note that p = -O.OOST GPa 
gives better temperature dependence of cja of tetragonal 
phase than p = 0.0 or p = —2.0 GPa, though cja is still 
slightly overestimated. 

We also calculated total-energy surfaces for zone-center 
distortions of PbTiOa and SrTiOa in Figs. [6] and [71 re- 
spectively. Note that SrTiOa is not a ferroelectric ma- 
terial and the polarizing zone-center distortion is not to 
be realized. However, these may be useful data for in- 
vestigating epitaxial constraint SrTiOa films where the 
polarization properties depend crucially on the lattice 
mismatch. In both cases, similar trends in energetics 
as those in BaTiOs can be seen: LDA results in shal- 
low double wells and GGA (PBE) results in rather deep 
wells, while GGA (Wu and Cohen) results in intermediate 
depths of the double well potentials. Again, LDA calcu- 
lations under negative pressures, —4.0 GPa for PbTiOs 
and —8.0 GPa for SrTiOs, give similar results to those 
of GGA (Wu and Cohen) . Equilibrium cubic lattice con- 
stants are compared in Tabled 
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FIG. 6. (Color onHne) (a)-(d) Total-energy surfaces for zone- 
center distortions of PbTiOs . (e) u-dependence of strain along 
[110] distortion, (f) "Direction" ^(n) of atomic displacements 
along [001] distortion. 



FIG. 7. (Color online) (a)-(d) Total-energy surfaces for zone- 
center distortions of SrTiOa. (e) u-dependence of strain along 
[110] distortion, (f) "Direction" ^(n) of atomic displacements 
along [001] distortion. 



TABLE III. Simulated cubic -f-)- tetragonal, tetragonal -H- or- 
thorhombic, and orthorhombic ^->- rhombohedral transition 
temperatures are compared to the previous MD simulations^ 
with LDA-based parameters and experimentally observed val- 
ues. Heating-up and cooling-down transition temperatures 
are averaged when corresponding transition has temperature 
hysteresis. 



XC functional and 


ortho. 


tetra. 


cubic 


effective negative pressure 


rhombo 


ortho. 


tetra. 


GGA (W&C), 0.0 GPa 


102 K 


160 K 


288 K 


GGA (W&C), -2.0 GPa 


117 K 


218 K 


408 K 


GGA (W&C), -O.OOST GPa 


103 K 


187 K 


411 K 


LDA, 0.0 GPa 


95 K 


110 K 


137 K 


LDA, -5.0 GPa 


210 K 


245 K 


320 K 


experiment (after Refs. 24 and 


3) 183 K 


278 K 


403 K 



IV. SUMMARY 

In this work, we have evaluated the improvement in 
the description of energy surface at K relevant to ferro- 
electricity in perovskite-based titanates with use of Wu- 
Cohen GGA functional for exchange correlation energy 
in first-principles density functional theoretical calcula- 
tions to correctly estimate lattice constants and crys- 
tallographic anisotropy such as c/a ratio. We have 
demonstrated that the new GGA (Wu and Cohen) 
functional based calculations and LDA calculations un- 
der certain negative pressures are capable of yielding 
fairly accurate and comparable total-energy surfaces of 
zone-center distortions for ABO3 perovskite-type ferro- 
electrics; BaTiOa, PbTiOs, and SrTiOs. We have shown 
that their polar structural distortions are highly sensi- 
tive to their lattice constants, and hence much of the 
improvement with Wu-Cohen functional comes from its 
ability to correctly estimate lattice constants. We note 
that the use of Wu-Cohen functional has little influence 
elastic parameters, the soft mode eigenvectors and mode 
effective charges which govern the long-range dipolar in- 
teractions. What is most affected are the terms of the 



cubic anisotropy, as reflected in the strain coupling term 
Biyy and anharmoinc terms, and hence in the relative 
energy well-depths of polar distortions along [001], [110], 
and [111] directions. This is not quite surprising as the 
motivation in development of the Wu-Cohen functional 
was to get structural properties such as c/a ratio and 
lattice constants with better accuracy. 

We then analyzed consequences of this improvement 
in energy functional to finite temperature ferroelectric 
transition, taking an example of BaTiOa. To this end, 
starting from calculations with GGA (Wu and Cohen) 
functional, we constructed a new parameter set for effec- 
tive Hamiltonian of BaTiOs employing the valley-tracing 
technique that effectively includes anharmonic coupling 
of the soft polar mode with higher energy polar modes. 
Comparing this and an LDA-based effective Hamiltonian 
with MD simulations, we find that the use of Wu-Cohen 
functional leads to a clear improvement in description of 
the highest temperature transition from cubic to tetrag- 
onal phase. We also confirmed that, as already men- 
tioned in Refs. [2a and [2^, the effect of thermal expan- 
sion, which is basically coming from the odd order en- 
ergy terms of atomic displacements and their coupling 
with strains, cannot be ignored as these materials ex- 
hibit strong electro-mechanical couplings. Secondly, ac- 
counting for thermal expansion approximately through 
the temperature-dependent effective negative pressure in 
effective Hamiltonian, a more realistic description of fer- 
roelectric phase transition can be obtained. 
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